EMA analyses

load packages

library(pacman)
pacman::p_load(tidyverse, brms, ggeffects, kableExtra, tidybayes, install = TRUE)

define aesthetics

palette = c("#3DC0E1", "#1985A1", "#4C5C68", "#00A896")
palette_group = c("#0C7BDC", "#D55E00")
palette_sd = c(palette[1], palette[3], palette[2])

plot_aes = theme_minimal() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        text = element_text(size = 16, family = "Futura Medium"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(color = "black"),
        axis.line = element_line(colour = "black"),
        axis.ticks.y = element_blank())

define functions

make_table = function(data) {
    data %>%
      broom.mixed::tidy(conf.int = TRUE) %>%
      filter(effect == "fixed") %>%
      mutate(term = gsub("\\(Intercept\\)", "intercept", term),
             term = gsub("regulation_expression", "signature expression (mindfulness)", term),
             term = gsub("active_weekon", "intervention week (on)", term),
             term = gsub("active_weekoff", "intervention week (off)", term),
             term = gsub("signal_count", "signal", term),
             term = gsub("SocialWeekendWeekend", "weekend", term),
             term = gsub("gender_f", "gender (female)", term),
             term = gsub(":", " x ", term),
             `b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high)) %>%
      select(term, `b [95% CI]`) %>%
      knitr::kable(digits = 2)
}

load and tidy data

• remove trials missing stimulus information
  • the following participants had technical errors during all mindfulness trials: sub-MURIP012, sub-MURIP063, sub-MURIP078
• remove trials in which the mean signal intensity of the brain map is +/- 3 SDs from the mean

task

task_data = read.csv("../data/task/cuereact_all062920.csv", stringsAsFactors = FALSE) %>%
  filter(grepl("rating", trial_type)) %>%
  select(pID, block_num, trial_num, trial_type, resp, rt, stim) %>%
  mutate(trial_cond = ifelse(grepl("friend3|friend4", trial_type), "regulation",
                      ifelse(grepl("friend1|friend2", trial_type), "up-regulation",
                      ifelse(grepl("mindful", trial_type), "regulation",
                      ifelse(grepl("nonalc", trial_type), "non-alcohol reactivity",
                      ifelse(grepl("rating_alc_react", trial_type), "reactivity", NA))))),
         trial_cond = factor(trial_cond, levels = c("non-alcohol reactivity", "reactivity", "regulation", "up-regulation"))) %>%
  rename("stimulus" = stim) %>%
  filter(!is.na(stimulus))

task_conditions = task_data %>%
  select(pID, trial_type) %>%
  unique() %>%
  filter(grepl("friend|mindful", trial_type) | pID %in% c("sub-MURIP012", "sub-MURIP063", "sub-MURIP078")) %>%
  mutate(condition = ifelse(grepl("friend", trial_type), "perspective", "mindfulness")) %>%
  select(-trial_type) %>%
  unique()

task_data_all = task_data %>%
  left_join(., task_conditions) %>%
  mutate(condition = ifelse(is.na(condition), "control", condition),
         trial_cond = factor(trial_cond, c("non-alcohol reactivity", "reactivity", "regulation", "up-regulation")))

task_mindfulness = task_data_all %>%
  filter(condition == "mindfulness")

dots

file_dir = "../data/dots"
file_pattern = "sub.*"
file_list = list.files(file_dir, pattern = file_pattern)

dots_tmp = data.frame()

for (file in file_list) {
  tmp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE, header = TRUE, sep = ",") %>%
                    extract(file, c("pID", "beta"), ".*(sub-MURIC[0-9]{3}|sub-MURIP[0-9]{3})/(beta_[0-9]{4}).nii") %>%
                    mutate(stimulus = as.character(stimulus)), error = function(e) message(file))
  
  dots_tmp = rbind(dots_tmp, tmp)
  rm(tmp)
}

dots = dots_tmp %>%
  left_join(., task_conditions) %>%
  mutate(condition = ifelse(is.na(condition), "control", condition)) %>%
  group_by(condition) %>%
  filter(!stimulus == "nan") %>%
  mutate(sd3_signal = 3 * sd(mean_signal, na.rm = TRUE),
         grand_mean_signal = mean(mean_signal, na.rm = TRUE),
         outlier = ifelse(mean_signal > grand_mean_signal + sd3_signal, 1,
                   ifelse(mean_signal < grand_mean_signal - sd3_signal, 1, 0))) %>%
  select(-sd3_signal, -grand_mean_signal) %>%
  filter(!outlier == 1)

classifier data

classifier_data = read.csv("../data/mvpa/logistic_stats.csv", stringsAsFactors = FALSE) %>%
  rename("pID" = subjectID) %>%
    mutate(predicted_factor = ifelse(predicted > .5, "regulate", "react"),
           predicted_factor = as.factor(predicted_factor),
           actual_factor = ifelse(actual == 1, "regulate", "react"),
           actual_factor = as.factor(actual_factor))

merge and prep for modeling

• filter out motion exclusions: sub-MURIC303
• filter out reactivity to non-alcoholic beverage trials

merged = task_data_all %>%
  full_join(., dots) %>%
  filter(!pID == "sub-MURIC303") %>%
  filter(trial_cond %in% c("reactivity", "regulation")) %>%
  mutate(trial_cond_recode = ifelse(trial_cond == "regulation", .5, -.5),
         trial_cond = as.factor(as.character(trial_cond)))

between = merged %>%
  select(pID, dot, trial_cond, trial_cond_recode, condition) %>%
  group_by(pID, trial_cond, trial_cond_recode, condition) %>%
  summarize(dot_between = mean(dot, na.rm = TRUE)) %>%
  group_by(trial_cond, condition) %>%
  mutate(dot_between_c = scale(dot_between, scale = FALSE, center = TRUE)) %>%
  group_by(trial_cond, condition) %>%
  mutate(sd_dot = sd(dot_between, na.rm = TRUE),
         dot_between_std = dot_between_c / sd_dot,
         dot_between_std_noc = dot_between / sd_dot) %>%
  select(-sd_dot)

reg_react_expression = between %>%
  select(pID, condition, trial_cond, dot_between_std) %>%
  unique() %>%
  mutate(trial_cond = sprintf("%s_expression", trial_cond)) %>%
  spread(trial_cond, dot_between_std)

within-person regulation expression only

ema = read.csv("../data/ema/cleaned_EMA.csv", stringsAsFactors = FALSE) %>%
  select(pID, groupID, Condition, contains("signal"), Session.Name, active_week, drinks_number, contains("Craving"), gender_f, NumberResponses, SocialWeekend, amount_self_c, freq_self_c, contains("React")) %>%
  mutate(craving_previous = lag(Craving_Alc),
         Condition = gsub("mindful", "mindfulness", Condition),
         pID = gsub("muric", "sub-MURIC", pID),
         pID = gsub("muri", "sub-MURIP", pID),
         active_week = gsub("active", "on", active_week),
         active_week = gsub("control", "off", active_week),
         time_of_day = ifelse(grepl("Morning", Session.Name), "morning",
                       ifelse(grepl("Evening", Session.Name), "evening", NA))) %>%
  rename("condition" = Condition) 

ema_within_intervetion = ema %>%
  left_join(., reg_react_expression) %>%
  filter(!condition == "control") %>%
  group_by(pID) %>%
  mutate(drinks_number = scale(drinks_number, center = TRUE, scale = FALSE),
         Alc_React_Mindful = scale(Alc_React_Mindful, center = TRUE, scale = FALSE),
         Alc_React_Perspective = scale(Alc_React_Perspective, center = TRUE, scale = FALSE),
         mindful_scale = scale(Alc_React_Mindful, center = TRUE, scale = TRUE),
         perspective_scale = scale(Alc_React_Perspective, center = TRUE, scale = TRUE),
         response_scale = ifelse(condition == "mindfulness", mindful_scale, perspective_scale),
         craving_previous = scale(craving_previous, center = TRUE, scale = TRUE),
         craving = scale(Craving_Alc, center = TRUE, scale = TRUE))

ema_within = ema_within_intervetion %>%
  filter(condition == "mindfulness")

ema_between = ema %>%
  group_by(pID, condition) %>%
  summarize(drinks_between = mean(drinks_number, na.rm = TRUE),
            responses_between = mean(Alc_React_Mindful, na.rm = TRUE),
            craving_between = mean(Craving_Alc, na.rm = TRUE)) %>%
  filter(condition == "mindfulness") %>%
  left_join(., reg_react_expression)

visualize

intervention week effects

ema_within %>%
  select(pID, drinks_number, mindful_scale, craving, active_week) %>%
  gather(key, value, -pID, -active_week) %>%
  ggplot(aes(active_week, value, color = key)) +
  stat_summary(aes(group = interaction(pID, key)), fun = "mean", geom = "line", color = "grey", size = .2) +
  stat_summary(aes(group = key), fun = "mean", geom = "line") +
  stat_summary(fun.data = "mean_cl_boot") +
  facet_grid(~key) +
  plot_aes +
  theme(legend.position = "none")

correlations

mindful responses & drinking

Association between how mindful you were when you encountered alcohol and how much you drank (morning signal = reported evening drinking, evening signal = reporting morning/afternoon drinking)

ema_within %>%
  ggplot(aes(mindful_scale, drinks_number, color = active_week, fill = active_week)) +
  geom_point(size = .5, alpha = .4) +
  geom_smooth(aes(group = interaction(pID, active_week)), method = "lm", se = FALSE, size = .2) +
  geom_smooth(method = "lm", alpha = .2) +
  scale_color_manual(values = palette) +
  scale_fill_manual(values = palette) +
  plot_aes

craving (lagged) & drinking

E.g. association between morning craving and morning/afternoon drinking reported in the evening

ema_within %>%
  ggplot(aes(craving_previous, drinks_number, color = active_week, fill = active_week)) +
  geom_point(size = .5, alpha = .4) +
  geom_smooth(aes(group = interaction(pID, active_week)), method = "lm", se = FALSE, size = .2) +
  geom_smooth(method = "lm", alpha = .2) +
  scale_color_manual(values = palette) +
  scale_fill_manual(values = palette) +
  plot_aes

craving (lagged) & responses

E.g. association between morning craving and morning/afternoon mindful responses to alcohol reported in the evening

ema_within %>%
  ggplot(aes(craving_previous, mindful_scale, color = active_week, fill = active_week)) +
  geom_point(size = .5, alpha = .4) +
  geom_smooth(aes(group = interaction(pID, active_week)), method = "lm", se = FALSE, size = .2) +
  geom_smooth(method = "lm", alpha = .2) +
  scale_color_manual(values = palette) +
  scale_fill_manual(values = palette) +
  plot_aes

responses & craving

E.g. association between morning/afternoon mindful responses to alcohol reported in the evening and craving reported in the evening

ema_within %>%
  ggplot(aes(mindful_scale, craving, color = active_week, fill = active_week)) +
  geom_point(size = .5, alpha = .4) +
  geom_smooth(aes(group = interaction(pID, active_week)), method = "lm", se = FALSE, size = .2) +
  geom_smooth(method = "lm", alpha = .2) +
  scale_color_manual(values = palette) +
  scale_fill_manual(values = palette) +
  plot_aes

brms models

drinking ~ active_week * regulation_expression

plot

regression

vals = seq(-2, 2, .2)
predicted = ggeffects::ggpredict(drinking_week, c("regulation_expression [vals]",  "active_week")) %>%
  data.frame()

predicted %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2, color = NA) +
  geom_line(size = 1) +
  scale_color_manual(name = "intervention week", values = palette[3:4]) +
  scale_fill_manual(name = "intervention week", values = palette[3:4]) +
  labs(y = "within-person number of drinks\n", x = "\nsignature expression on mindfulness trials (SD)") +
  plot_aes + 
  theme(legend.position = "top")

point range

predicted = ggeffects::ggpredict(drinking_week, c("regulation_expression [0, 1]", "active_week")) %>%
  data.frame() %>%
   mutate(x = ifelse(x == 1, "+1 SD", "mean"),
          x = factor(x, levels = c("mean", "+1 SD")),
          group = recode(group, "on" = "active", "off" = "control"),
          group = factor(group, levels = c("control", "active")))

predicted %>%
  filter(x == "mean") %>%
  ggplot(aes(group, predicted, color = x, group = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(0.05), size = 1.5) +
  geom_line(position = position_dodge(0.05), size = 1.5) +
  scale_color_manual(name = "mindfulness signature expression", values = palette) +
  scale_x_discrete(expand = c(.1, .1)) +
  coord_cartesian(ylim = c(-.3, .3)) +
  labs(y = "within-person number of drinks\n", x = "\nintervention week") +
  plot_aes +
  theme(legend.position = "top")

predicted %>%
  ggplot(aes(group, predicted, color = x, group = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(0.05), size = 1.5) +
  geom_line(position = position_dodge(0.05), size = 1.5) +
  scale_color_manual(name = "mindfulness signature expression", values = palette) +
  scale_x_discrete(expand = c(.1, .1)) +
  coord_cartesian(ylim = c(-.3, .3)) +
  labs(y = "within-person number of drinks\n", x = "\nintervention week") +
  plot_aes +
  theme(legend.position = "top")

summary

make_table(drinking_week)

term	b [95% CI]
intercept	0.20 [0.06, 0.34]
intervention week (on)	-0.07 [-0.22, 0.08]
signature expression (mindfulness)	0.11 [0.01, 0.21]
signal	-0.01 [-0.01, -0.00]
intervention week (on) x signature expression (mindfulness)	-0.23 [-0.38, -0.08]

simple slopes

hypothesis(drinking_week,
           c(expression_mean = "active_weekon = 0",
             expression_1sd = "active_weekon + active_weekon:regulation_expression = 0"))$hypothesis %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", Estimate, CI.Lower, CI.Upper)) %>%
  select(Hypothesis, `b [95% CI]`) %>%
  knitr::kable()

Hypothesis	b [95% CI]
expression_mean	-0.07 [-0.22, 0.08]
expression_1sd	-0.30 [-0.50, -0.09]

responses ~ active_week * regulation_expression

plot

regression

vals = seq(-2, 2, .2)
predicted = ggeffects::ggpredict(response_week, c("regulation_expression [vals]",  "active_week")) %>%
  data.frame()

predicted %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2, color = NA) +
  geom_line(size = 1) +
  scale_color_manual(name = "intervention week", values = palette[3:4]) +
  scale_fill_manual(name = "intervention week", values = palette[3:4]) +
  labs(y = "within-person mindful response rating (SD)\n", x = "\nsignature expression on mindfulness trials (SD)") +
  plot_aes +
  theme(legend.position = "top")

point range

predicted = ggeffects::ggpredict(response_week, c("regulation_expression [0, 1]", "active_week")) %>%
  data.frame() %>%
   mutate(x = ifelse(x == 1, "+1 SD", "mean"),
          x = factor(x, levels = c("mean", "+1 SD")),
          group = recode(group, "on" = "active", "off" = "control"),
          group = factor(group, levels = c("control", "active")))

predicted %>%
  filter(x == "mean") %>%
  ggplot(aes(group, predicted, color = x, group = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(0.05), size = 1.5) +
  geom_line(position = position_dodge(0.05), size = 1.5) +
  scale_color_manual(name = "mindfulness signature expression", values = palette) +
  scale_x_discrete(expand = c(.1, .1)) +
  coord_cartesian(ylim = c(-.5, .65)) +
  labs(y = "within-person mindful response rating (SD)\n", x = "\nintervention week") +
  plot_aes +
  theme(legend.position = "top")

predicted %>%
  ggplot(aes(group, predicted, color = x, group = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(0.05), size = 1.5) +
  geom_line(position = position_dodge(0.05), size = 1.5) +
  scale_color_manual(name = "mindfulness signature expression", values = palette) +
  scale_x_discrete(expand = c(.1, .1)) +
  coord_cartesian(ylim = c(-.5, .65)) +
  labs(y = "within-person mindful response rating (SD)\n", x = "\nintervention week") +
  plot_aes +
  theme(legend.position = "top")

summary

make_table(response_week)

term	b [95% CI]
intercept	-0.49 [-0.73, -0.27]
intervention week (on)	0.44 [0.19, 0.69]
signature expression (mindfulness)	-0.10 [-0.24, 0.06]
signal	0.01 [0.00, 0.02]
intervention week (on) x signature expression (mindfulness)	0.26 [-0.01, 0.51]

simple slopes

hypothesis(response_week,
           c(expression_mean = "active_weekon = 0",
             expression_1sd = "active_weekon + active_weekon:regulation_expression = 0"))$hypothesis %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", Estimate, CI.Lower, CI.Upper)) %>%
  select(Hypothesis, `b [95% CI]`) %>%
  knitr::kable()

Hypothesis	b [95% CI]
expression_mean	0.44 [0.19, 0.69]
expression_1sd	0.71 [0.36, 1.03]

mediation

hypothesis(
  fit_brm,
  'b_drinksnumber_mindful_scale * b_mindfulscale_active_weekon + cor_pID__mindfulscale_active_weekon__drinksnumber_active_weekon * sd_pID__mindfulscale_active_weekon * sd_pID__drinksnumber_mindful_scale = 0',
  class = NULL,
  seed  =  6523
)

## Hypothesis Tests for class :
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (b_drinksnumber_m... = 0    -0.31      0.14    -0.59    -0.03         NA
##   Post.Prob Star
## 1        NA    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.